Ab initio and periodic DFT investigation of hydrogen storage on light metal- decorated MOF-5

The effect of light metal (M = Li, Be, Mg, and Al) decoration on the stability of metal organic framework MOF-5 and its hydrogen adsorption is investigated by ab initio and periodic density functional theory (DFT) calculations by employing models of the form BDC:M2:nH2 and MOF-5:M2:nH2, where BDC stands for the benzenedicarboxylate organic linker and MOF5 represents the primitive unit cell. The suitability of the periodic DFT method employing the GGA-PBE functional is tested against MP2/6-311+G* and MP2/cc-pVTZ molecular calculations. A correlation between the charge transfer and interaction energies is revealed. The metal-MOF-5 interactions are analyzed using the frontier molecular orbital approach. Difference charge density plots show that H2 molecules get polarized due to the charge generated on the metal atom adsorbed over the BDC linker, resulting in electrostatic guest-host interactions. Our solid state results show that amongst the four metal atoms, Mg and Be decoration does not stabilize the MOF-5 to any significant extent. Li and Al decoration strengthened the H2-MOF-5 interactions relative to the pure MOF-5 exhibited by the enhanced binding energies. The hydrogen binding energies for the Liand Al-decorated MOF-5 were found to be sensible for allowing reversible hydrogen storage at ambient temperatures. A high hydrogen uptake of 4.3 wt.% and 3.9 wt.% is also predicted for the Liand Al-decorated MOF-5, respectively. Depletion of fossil fuels as well as environmental concerns has led to a widespread research for finding new and renewable energy sources. Hydrogen (H2), with its high chemical energy content (142 MJ/kg), that is at least three times larger than that of other chemical fuels (for example, the equivalent value for liquid hydrocarbons is 47 MJ/kg), abundance and nonpolluting nature has generated tremendous interest [1]. Hydrogen economy offers an alternative to our dependence on fossil fuels and provides an environmentally clean source of energy. However, the application of hydrogen economy for practical applications requires efficient H2 storage with high gravimetric and volumetric density. In comparison to storage of H2 in a high-pressure tank and cryogenic liquefaction of hydrogen, hydrogen adsorption in solid materials is much safer [2]. A number of research efforts have been made to find new storage materials which are capable of exhibiting reversible H2 storage at high densities, reasonably fast adsorption/desorption kinetics and long life time [3]. However, till date no material has been found to have all the above stated properties. Solid state hydrogen adsorption is currently pursued by two main strategies: strong chemical associations (chemisorption) and weak dispersive interactions (physisorption). As an example, the chemisorption approach applies to the class of simple and complex metal hydrides. Although these materials show high H2 storage, they suffer from large dehydrogenation activation barriers. Using periodic density functional theory (DFT) calculations, Kelkar et al. [4] have shown that doping different phases of MgH2 with Al and Si can result in a lowering of its activation energy barriers associated with the direct desorption of H2 from its (001) surface. With regard to physisorption of H2, materials with large surface areas and low densities, such as metal organic frameworks (MOFs), are attractive. MOFs are crystalline hybrid organic/inorganic nanoporous materials which can adsorb guest molecules, and thus are potential candidates for H2 storage [5]. There are predominantly Van der Waals interactions between the physisorbed H2 molecules and the host MOFs. Thus, the resultingH2 adsorption enthalpies are typically in the range of 2.26e5.2 kJ/mol [6-8]. As a result, MOFs exhibit fast adsorption/desorption kinetics. The drawback is that significant hydrogen adsorption can only be achieved at low temperatures ofw77 K. A study by Lochan and Head-Gordon [9] has shown that the ideal H2 binding energy range is 28-40 kJ/ mol at room temperature andw30-42 kJ/mol ~50 C. Thus, clearly there is a significant need to enhance the H2 binding energy (Ebinding) as well as increase H2 adsorption and storage capacity of MOFs in order to render them useful for practical and mobile applications. The various approaches to improve hydrogen adsorption in MOFs include the optimization of pore size [10], adding another guest molecule in the pore of crystal [11] and catenation [10]. Li et al. [12] have recently exploited the idea of optimization of the pore size. They have shown that MOF-5 synthesized by slovthermal approach has higher surface area and pore volume, than by the conventional methods of direct mixing or slow diffusion of triethyamine. To enhance the Ebinding, introduction of electron donating groups and inclusion of open metals sites [13] have been attempted. It is necessary to develop an understanding of the fundamental interactions involved in hydrogen adsorption in order to be able to tune H2 storage in MOFs. Ab inito and DFT studies on the molecular linker alone have been found useful to predict the hydrogen binding sites and the orientation and to determine the binding energies. The first theoretical study on MOF-5 was carried out by Hubner et al. [14] by investigating simplistic model systems, namely, C6H6X (X= H, F, OH, NH2, CH3 and CN), naphathalene, anthracene, coranene, terethalicacid and dilithiumterephthalete. Their calculations carried out by the second-order Møller-Plesset (MP2) method revealed that the Ebinding increased with increase in ring size from 3.91 kJ/mol for benzene to 4.28 kJ/mol for naphthalene. In a similar work Sagara et al.[15] studied both the organic linker part and zinc oxide part of MOF5 using MP2 method. They concluded that H2 preferably binds between two zinc clusters. In another study, Buda et al. [16] have illustrated that the top site had a higher Ebinding than the edge site of BDC linker in MOF-5. In a series of landmark papers from 2007 to 2009, Han et al.[17-19] showed that accurate force fields can be fitted to such MP2 or CCSD(T) calculations and then used in Grand Canonical Monte Carlo (GCMC) calculations to determine the loading of H2 as a function of pressure and temperature. They demonstrated that this very accurately reproduces the experimental loading curve for MOF-5 and other MOF and COF systems. Mulder et al. [20] and Mueller and Ceder [21] carried out a full periodic DFT study of MOF-5 and found five symmetrically unique adsorption sites, their results indicated that the sites with the strongest interaction with hydrogen are located near the Zn4O clusters. The recent advances in studying MOFs and covalent organic frameworks (COFs) using quantum calculations, MonteCarlo (MC) simulations and molecular dynamics simulations have been critically reviewed in Ref. [22]. The authors have also discussed the different strategies being pursued for improving hydrogen storage in these materials. One such interesting method involves decoration of the organic linker part of the MOFs with a metal. Previously many theoretical studies have been focused on the interaction of H2 molecules with transition metal ions [2326]. Noteworthy is the investigation by Gagliardi and Pyykko [27] which showed that bare neutral and cationic transition metals such as Cr, Mn, W, Mo, V and Ti can bind with as many as six H2 molecules. Maark et al. [28] have studied the effect of metal ions, Na, Mg, Al and Be using benzene as a simple model of MOF-5. A well-known study by Chandrakumar and Ghosh [29] showed that s-block metal cations could bind with even more number of H2 molecules and form complexes of the type M(H2)8 complexes with high binding energies. The binding energies for these light metal cations were calculated to be between 56.484 kJ/mol to 753.1 kJ/mol. However, manufacturing such bare metal cations in reality is difficult. Alternatively the metal atoms can acquire charge via charge transfer by adsorption on surfaces with high electron affinity. Recently, activation of nanomaterials has also been attempted via spillover effect. Research studies have exhibited that at moderate temperature and pressure conditions, doping of nanomaterials like nanotubes [30], nanofibers [31] and nanospheres [32] increases the Ebinding of these materials. But transition metal based materials show very high binding energy [33-37]. Zhao et al. [33], through their studies on organometallic buckyballs, found that the calculated Ebinding for transitions metals was w0.3 eV. In another well known study of transition metal decoration by Yildirim et al. [34] illustrated that Ti-decorated nanotubes can adsorb upto 8 H2 per Ti but with a high binding energy of 4.3 eV/H2. Due to these high Ebinding, these materials cannot be used for reversible H2 storage. Furthermore, transition metals have a tendency to bind among each other thereby decreasing the H2 storage capacity [38]. This has been found to be true in case of embedding in MOFs [39-44] or attachments by strong covalent bonds with host materials [37,45]. A possible solution of this problem is to use elements with low cohesive energy, like s-block metals. The light metal decoration approach has been investigated theoretically on several systems including carbon nanotubes (CNTs), graphene and fullerenes [33,35,46-52]. In a particularly interesting work, Sun et al. [48], using density functionaltheory, illustrated that Li-coated fullerenes do not suffer fromthe clustering problem of transition metal based materials discussed above. They found that an isolated Li12C60 cluster where Li atoms are capped onto the pentagonal faces of the fullerene was not only very stable but could also store up to 120 hydrogen atoms in molecular form with a binding energy of 0.075 eV/H2. A similar strategy is now being widely considered for MOFs. Blomqvist et al. [53] have studied Li-decoratedMOF-5 by ab initio periodic DFT calculations and shown that each Li adsorbed over each benzene ring in MOF-5 can cluster up to three H2 molecules with binding energies between 12.0-18.0 kJ/mol H2. It was predicted through ab initio molecular dynamics simulations that the system would exhibit a hydrogen uptake of 2.9 wt.% and 2.0 wt.% at the high temperatures of 200 K and 300 K, respectively. Han and Goddard III [54] have predicted gravimetric H2 adsorption isotherms for several pure and Li-doped MOFs at 300 K using an ab initio based GCMC simulation. They found that Li-doped MOF-5 showed significant H2 binding at 300 K and 100 bar. Another interesting result of their calculations was the high H2 uptake of 6.47 wt.% predicted for MOF-C30 at these conditions. Recently, Venkataramanan et al. [55] have explored isoreticular MOFs with different metals M=( Fe, Cu, Co, Ni and Zn) through DFT and predicted that Li doping will be possible only in Zn-based MOFs. Ab initio calculations have demonstrated that Ca can decorate organic linkers of MOF-5 with a binding energy of 1.25 eV and that the resulting hydrogen storage can be as high as 4.6 wt.% [56] Furthermore, the H2 binding was found to be significantly stronger than the Van der Waals interactions. In a current study Zou et al. [57] have studied hydrogen storage in Ca-decorated, B-substituted MOFs through first-principles electronic structure calculations. Substitution of B atoms in place of C atoms in the benzene ring of the BDC linker increased strength of interaction of Ca with the linker as well as the interaction of H2 molecule with Ca. A hydrogen uptake of 8 H2 per Ca-decorated linker was found with binding energy of 20 kJ/mol rendering it suitable for reversible H2 storage under ambient conditions. Many attempts have been made to determine the major contributing factor, the London dispersion interactions (LDI) versus interactions due to the electrostatic potential of the host materials, for H2 adsorption in pure MOFs. Bordiga et al. [7] applied a combination of IR spectroscopy and ab initio calculations to exhibit that the adsorptive properties of MOF-5 are mainly due to the (i) dispersive interactions with the internal wall structure and (ii) weak electrostatic forces associated with O13Zn4 clusters. Kuc et al. [58,59] modeled IRMOF-1 by individually examining its cornerposts andzlinkers. Their calculations revealed that the physisorption of H2 in MOFs is mainly due to LDI between linkers and connectors with hydrogen while the host-guest induced electrostatic interactions were unimportant as the charge separation in the MOF is not large enough to induce significant dipole moments in H2. It is further concluded that it is necessary to use correlation methods and large basis sets for correctly describing the LDI in these systems. Therefore, methods based on DFT which are unable to correctly describe dispersion forces and Van der Waals interactions face challenges in studying pure MOFs. But it is also known that ab initio methods such as MP2 and CCSD(T) are not practical for application to the full structure. Motivated by the above facts, herein we have studied metal (M)-decorated MOF-5, where M= Li, Be, Mg and Al, using both periodic DFT based calculations employing the PBE functional and molecular quantum chemical computations with MP2 methods to test the accuracy of the former for predicting the structures and H2 binding energies of the metal-decorated MOF-5s. Further, we have analyzed which light metal atom will be the best for increasing the strength of H2 adsorption and amount of H2 storage in MOF-5 for making it suitable for ambient H2 storage applications. To this end we have examined the optimized geometries, atomic charges, interaction energy of the metal atoms with MOF-5 and the hydrogen binding energies. 2. Computational methodology It has been reported by Mueller and Ceder [17] that studies on a truncated part of MOF like BDC linker alone can yield misleading results. In order to investigate this effect, in this work we have performed calculations on two sets of systems: (i) the full primitive cell of MOF-5 and (ii) the molecular benzenedicarboxylate(BDC) linker. In the paper the pure and metal decorated MOF-5 unit cells are referred to as MOF-5 and MOF5:M2, respectively and the corresponding molecular counterparts are referred as BDC and BDC:M2, respectively, where M= Li, Be, Mg, and Al. The interaction energy (ΔE ) of the metal atom with MOF5 is calculated as follows: E(MOF Δ -5:M) =1/2[Etot (MOF-5:M2) – 2*Etot(M) – Etot(MOF-5)] ..........(1) Different numbers of H2 molecules were allowed to interact with Li-decorated BDC and M-decorated MOF-5. These are represented as BDC:Li2:nH2 and MOF5:M2:nH2, respectively, where n is the number of H2 molecules adsorbed. The corresponding H2 binding energies per H2 are calculated by the equations: Ebinding(BDC:M2:nH2) =1/n[E(BDC:M2:nH2) − E(BDC:M2 ) – n E(H2 )] ..........(2) Ebinding(MOF-5:M2:nH2) =1/n[E(MOF-5:M2:nH2) − E(MOF-5:M2 ) – n E(H2 )] ..........(3) Density functional theory (DFT) [60] is capable of scaling well with system size. It, however, unfortunately fails to calculate magnitude of weak interactions accurately because these interactions depend strongly on electron correlation. Within the local density approximation (LDA) these weak interactions get overestimated [61-64]. In comparison the generalized gradient approximation within the Perdew, Bruke and Ernzerhof (PBE) [65] exchange correlation functional is known to give reasonably accurate results [61,63]. However, it is still important to test the reliability of using DFT for the systems under study. For this purpose we have implemented three computational methods: gas-phase and periodic DFT calculations with the PBE functional using a plane wavepseudopotential approach and molecular calculations with MP2 and other DFT functionals and 6-311+G* and cc-pVTZ basis-sets. All the DFT planewave-pseudopotential based calculations have been performed using the Vienna ab initio simulations package (VASP) [66,67]. We have used the projector augmented wave (PAW) [68,69] approach to evaluate all the properties. An energy cut-off of 520 meV was employed throughout our periodic DFT calculations. All geometry optimizations were carried out without any geometry constraint. All forces were calculated using HellmannFeynman theorem. Geometries considered optimized with maximum force found smaller than 0.01 eV/Ǻ. Real space projections were used to evaluate PAW character of wave-functions. The valence states for all the potentials used here are 4s and 3d for Zn, 2s and 2p for both O and C and 1s for H. The primitive cell of MOF-5 was built from the crystal structure of MOF-5 that was taken from experimental data [70] and a 2 2 2 k-point grid generated by Monkhorste Pack scheme [71] was applied to it. For the gas phase calculations of pure and decorated BDCs using VASP a large cell size of 20 20 20 was chosen to avoid interaction between successive BDC molecules. These calculations were carried out at gamma-point. For obtaining the atomic charges Bader analysis [72] was performed. Our quantum chemical calculations with MP2 and other DFT functionals with the 6-311+G* and cc-pvtz basis sets were carried out with the Gaussian 09 package [73]. For all super-molecules basis set superposition error (BSSE) correction was evaluated by the counterpoise method [74]. The atomic charges were derived by the CHELPG ESP method [75]. All our calculations are 0 K calculations and do not include zero-point energy contribution. 3. Results and discussion 3.1. Structure Herein we begin with a comparison of the structures of BDC and BDC:M2 (M = Li, Al, Be and Mg) systems optimized from MP2/6-311+G* and GGA-PBE/plane wave methods. The MP2 optimized geometry of the model BDC linker is displayed in Fig. 1. Selected bond distances of the bare and metal decorated organic linker calculated from the two methods are presented in Table 1. It can be seen that in all cases there is a reasonable agreement between the structural parameters obtained from the two methods. These results suggest that using the PBE exchange correlation functional is reliable for structural investigations of this class of molecules. The effect of adsorption of metal on BDC can be understood by comparing the corresponding bond distances in Table 1. Clearly, except for C2-C3 bond distances predicted by MP2 method for BDC:Li2 and BDC:Al2, there is a marked increase in all other C1-C2 bond lengths in the different systems obtained from both the GGA-PBE and MP2 methods. Overall the C2-C3 elongation is most pronounced in case of BDC:Be2. The distance of metal atom from center of mass of BDC (BDCcom) is also listed in Table 1 and it exhibits the order: BDC:Be2 < BDC:Li2 < BDC:Mg2 < BDC:Al2. To determine the effect of going from a gas phase to a solid-state calculation on the structural parameters ofMOF5 we performed periodic DFT calculations with the PBE-GGA functional on the full primitive unit cell of MOF-5 consisting of 106 atoms. The optimized structure thus obtained is shown in Fig. 2. Table 2 lists the selected structural parameters calculated from this work, a previous theoretical study [21] and experiment [5]. It also presents the ratio of our calculated values with respect to the experimental results for each structural parameter of MOF-5. Our results show strong agreement with both theory and experiment. The latter is depicted by the high average ratio of 0.995. In the table we have also given the average C-C bond lengths in the BDC linker obtained from our GGA-PBE/plane wave and MP2/6-311+ G* calculations. A comparison with the experimental data shows that these values do not show as good an agreement. This tells us that for an accurate study of the structures of MOF-5 the primitive unit cell GGA-PBE computations are most suited. As a next step light metals were decorated on the top of the benzene ring (on both its sides) of MOF-5. The optimized primitive unit cells of the Li, Mg, Be and Al decorated-MOF-5s are showcased in Fig. 3. Following metal decoration different numbers of H2molecules were allowed to interact with MOF5. Selected structural parameters of MOF-5:M2, MOF5:Li2:nH2 and MOF-5:Al2:nH2 are listed in Table 3. The labeling of C atoms in the organic linker part of MOF-5 is the same as shown for BDC in Fig. 1. The table shows that the C-C bond lengths seem to remain unaffected in case of Mg. But in all other MOF-5:M2 systems, there is an elongation of C2-C3 bond and a contraction of the C1-C2 and C3-C3 bonds relative to the pure MOF-5. A look at the distances of the metal atoms to the COM of the BDC in the decorated MOF-5s reveals that Mg is more than twice the distance away from the BDC than the other atoms. The large values of Mg-BDCcom (3.754A° ) suggest that there is little interaction between Mg and MOF-5. It is further evident from Table 3 that addition of H2 molecules does not change the MOF-5 structure much. But as number of H2 molecules increases the M-BDCcom distance and the highest H2eM distance also increase along with a slight decrease in HeH bond length in the adsorbed H2 molecules suggesting a weakening of interactions with MOF-5. Though the M-BDCcom distances in MOF-5:Li2:nH2 and MOF5:Al2:nH2 for the same n are similar, the H2 molecules are further away from Al than Li implying a relatively weaker interaction with Al. On this basis we expect the Al-decorated MOF-5 to exhibit a lesser H2 binding energy than Li-decorated MOF-5. Comparing Table 3 to the values for BDC:M2 computed using the MP2/6-311+G* method presented in Table 1, it can be seen that the C-C bond lengths show a wide variation. For all M the C1-C2 and C3-C3 bonds in MOF-5:M2 are shorter and the C2-C3 bonds relatively longer than predicted for the BDC:M2 molecular systems. The magnitudes M-BDCCOM also differ for the solid state and gas phase systems. The MBDCCOM order in MOF-5:M2 is Be < Al Al > Mg > Be while the exhibited order is Li > Al > Be > Mg in MOF-5:M2 and Li > Mg > Al >Be in BDC:M2. Recalling the discussion in Section 3.1, adsorption of Mg on MOF-5 had little effect on its structural parameters. The Mg-BDCcom distance in MOF-5:Mg2 was w70% larger than in BDC:Mg2. Clearly, periodic GGA-PBE calculations predict a nonionic interaction between Mg and MOF-5 and so a direct comparison of charge transfer in the two systems would be incorrect. Removing Mg from the above series both calculation methods yield the same order of Q(M). Baker and Head-Gordon [76] have studied the interaction of Li with other polycyclic aromatic hydrocarbons with various electronic structure methods and observed an artificial electron transfer in some cases when using DFT due to the selfinteraction error. They found that for the ionic system Li(C6H6) the B3LYP and BLYP density functionals produced partial charge transfer of only 0.69 and 0.46 respectively for as opposed to 0.88 with the CCSD method. However, in our analysis the atomic charge for Li in MOF-5:Li2 derived from MP2 and GGA-PBE methods illustrated good agreement. For the other metals Table 4 shows that compared to MP2, DFT overestimates the charge transfer in MOF-5:Be2 and MOF-5:Al2 by nearly 50% and underestimates by the same percentage for MOF-5:Mg2. Furthermore, Bader charges have previously been shown to be substantially in error in a study of charge quadrupole interactions for H2 adsorption and diffusion in CuBTC [77]. Thus, whether the interaction energies between the metal atoms and MOF-5 conform to the amount of charge transfer predicted needs to be explored. The interaction energies (ΔE(MOF-5:M)) between the metal atom and MOF-5 in MOF-5:M2 calculated using Eq. (1) are also presented in Table 4. The calculated values of MOF-5:Be2 and MOF-5:Mg2 are much smaller in magnitude than those of MOF-5:Li2 and MOF-5:Al2. The ΔEs of MOF-5:Mg2, MOF-5:Li2 and MOF-5:Al2 are in accordance to the GGA-PBE/plane-wave calculated atomic charges, while ΔE of MOF-5:Be2 can only be a result of low charge transfer as predicted by MP2/6-311G* for BDC:Be2. This discussion suggests that except for Be decorated MOF-5, the GGA-PBE functional within the planewavepseudopotential approach can be applied to accurately model the charge transfer between the metal atoms Li, Mg and Al and MOF-5. A negative ΔE implies that the interaction between the metal atom and MOF-5 is exothermic and that the decoration stabilizes the MOF5. A positive value tells that interaction is endothermic and that metal decoration destabilizes the system. All the ΔEs in Table 4 are negative but those of MOF5:Be2 and MOF-5:Mg2 are nearly close to zero while MOF-5:Li2 and MOF-5:Al2 are considerably larger. Therefore, at room temperature Be and Mg may not remain adsorbed to MOF-5. These results suggest that out of the four metals investigated herein only Li and Al will stabilize MOF-5 to a significant extent and hence, H2 interactions are investigated only in the Liand Al-decorated systems. 3.3. Nature of metal-MOF-5 and metal-H2 interactions In this section we first correlate the nature of interactions between the metal atoms and MOF-5 to the observed structural parameters and interaction energies. It is understood that the effective strength of the metal-MOF-5 interaction will be a net result of the competing attractive metal-p-arene type interaction (due to the charge generated on the metal) and the mutual repulsion between the total electron densities of the metal atoms and benzene rings of MOF-5. To gain insight into the attractive interactions we have used the Frontier Molecular Orbital approach. It is known that the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) are the best global descriptors of reactivity [78,79]. In Fig. 4(a)e(e) we have illustrated the HOMOs of the four metal atoms and the LUMO of pure MOF-5. The HOMOs of the resulting metal decorated MOF-5s are displayed in Fig. 5(a)e(d). It can be seen that the HOMO of Li, Be and Mg is s type and of Al is p type. In MOF-5 the doubly degenerate unoccupied p-type orbital of benzene splits into two components because of the two carboxylic groups attached to it. The orbital component with a large density on the two ipso carbon atoms interacts constructively with the p* orbitals of the carboxylic groups, forming the MOF-5 LUMO (Fig. 4(e)) which is not completely delocalized over the benzene ring due to nodal planes. When the metal atoms are adsorbed over the BDC linker the HOMOof the metal atoms will interact with the LUMO of MOF-5. The HOMOs of the resulting metal-decorated MOF-5s are shown in Fig. 5(a)e(d). It can be seen that in case of MOF-5:Be2 and MOF-5:Mg2 HOMO the electron density is localized around the metal atoms showcasing that they have little interaction with MOF-5. So at room temperature these metal atoms may get desorbed from the metal-organic framework. In comparison the Liand Al-decorated MOF-5 HOMO are of the p type, delocalized around the benzene rings, similar to the LUMO of MOF-5 illustrating a transfer of electron density from the metal to the organic linker. Further, as there is no directional electron density between the metal and BDC, it can be inferred that the metal-MOF-5 interaction is ionic in nature and of the form M + δ (BDC)δwhere + δ is the charge generated on the metal due to electron transfer. The greater cloud distribution around the carbon atoms in the Lidecorated MOF-5 than in the Al-decorated MOF-5 is indicative of a greater electron transfer and hence, a stronger interaction of MOF-5 with Li than with Al. Based on the above discussion, the strength of the metal-parene attraction in MOF-5:M2 systems decreases as follows: Li > Al >>Be > Mg. As the electron density depends on the number of electrons, the mutual repulsion between the electron densities of the metalandMOF-5will increase in the order: Li < Be < Mg < Al. Clearly, the ΔE of MOF-5:Li will be the largest. Though a greater repulsion is expected between electron densities of Al and MOF-5, the contribution of the attraction on account of charge transfer is much higher than Be and Mg and therefore, MOF-5:Al2 will exhibit the next highest ΔE. Both Be and Mg HOMOs have poor attraction for the MOF-5 LUMO, but Mg has larger number of electrons. Therefore, the MOF-5:Mg interaction will be the weakest. This conforms to the metalMOF-5 interaction energy order (MOF-5:Li< MOF 5:Al <<MOF-5:Be< MOF-5:Mg) predicted in Section 3.2. At this point we would like to recall that the atomic charge on Li and Al in MOF-5:Li2 and MOF-5:Al2 was predicted to be +0.90 and +0.66 a.u., respectively (see Table 4). Therefore, when H2 molecules are added to these decorated MOF-5s they interact not with a neutral metal atom but with a metal ion. The transfer of electrons from metal ions to the H2 molecules is energetically not favorable because removing electrons from a positive ion is difficult by virtue of its much higher second ionization potential compared to the first ionization potential. The H2 adsorption therefore, takes place by a charge polarization mechanism in which the metal ion polarizes the nearby hydrogen molecules. The bonding between the metal ion and H2 molecule thereby proceeds through a simple electrostatic interaction. For understanding the above mechanism of charge polarization, we derived charge density difference surfaces of MOF5:Li2:2H2 (MOF-5:Li2:4H2) relative to the charge densities of MOF-5:Li2 and two H2 molecules (four H2 molecules). The charge density calculation for MOF-5:Li2 and two H2 molecules (four H2 molecules) was carried out by keeping the constituent atoms at the same positions as in MOF-5:Li2:2H2 (MOF5:Li2:4H2) and using the unit cell dimensions of the hydrogenated decorated MOF-5. The charge density difference plots, thus obtained, are showcased in Fig. 6. In these plots, red regions imply electron loss and blue regions indicate electron gain. It is noticeable that the change in electron density is occurring only around the H atoms. Around each H atom there is a blue region toward the Li ion and a corresponding red region on the other site. In this manner this picture depicts that due to polarization by the Li ion there is a shift of electron density in the H2 causing part of the H atoms to become positively charged (red region) and another part to become negatively charged (blue region). The experimental bond length of H2 molecule is 0.741A°. The H-H distances in MOF 5:Li2:nH2 and MOF-5:Al2:nH2, as listed in Table 3, are between 0.750 and 0.756A° . Thus, in the Li and Al-decorated MOF-5 hydrogen molecules are bonded in molecular form. The slight increase in bond length is an artifact of the polarization of the H2 molecule by the metal ion. 3.4. Hydrogen binding energies and gravimetric storage In order to systematically investigate the suitability of DFT functionals for the calculation of H2 binding energies for the metal decorated MOF-5s we first compare the BSSE corrected Ebinding for the BDC:Li2:nH2 molecular models evaluated using MP2 and different DFT functionals available with Gaussian 09 and using the GGA-PBE functional with the periodic DFT code VASP. Ebinding calculations utilizing 6-311+G* are shown in Fig. 7(a). It can be seen from Fig. 7(a) that there is considerable difference in the Ebinding by MP2 method and by GGA-PBE (VASP). This figure also suggests that B3LYP and PBEPBE functional underestimates binding energy while PW91PW91, PBE1PBE and GGA-PBE (VASP) functionals overestimates the binding energy. To consider long range interactions 6-311+G* basis set does not have sufficient polarizations functions and leads to unphysical large basis set errors [80]. To avoid such unphysical large basis set error, we have used cc-pvtz basis set. The results for Ebinding are shown in Fig. 7(b). In comparison, the H2 binding energies computed using the GGA-PBE functional along with the plane wave approach show excellent agreement with ccpvtz/ MP2 results for each n number of H2 molecules adsorbed. B3LYP Functional is found to underestimate Ebinding, while PW91PW91, PBEPBE and PBE1PBE functionals are found to highly overestimate Ebinding. The Ebinding for PBEPBE and PW91PW91 functionals are found to mach up to second decimal digit. The optimized geometries of the corresponding n = 2-6 H2 adsorbed Li-decorated MOF-5s are shown in Fig. 8. The total H2 binding energy and binding energy per H2 molecule of MOF5:Li2:nH2 and MOF-5:Al2:nH2 from the solidstate calculations of the full primitive unit cells are listed in Table 5. Our Ebinding per H2 for MOF-5:Li2:nH2 (n = 2-6) show an excellent agreement with the previously reported theoretical results of Blomqvist et al. [53] also presented in Fig. 7(b), as well as with the MP2 values calculated for the BDC:Li2:nH2 models in Fig. 7(b). This agreement depicts that the solid-state GGA-PBE methodology can be used to calculate with good precision the Ebinding for the Liand Al-decorated MOF-5s studied in this work. It is evident from both Table 5 and Fig. 7(b) that Ebinding decreases as the number of hydrogen molecules increases. This decrease is a result of the increased repulsion between the H2 molecules, the reduced sigma donation of the H2 molecules, and the consequent reduced effective charge on the metal atoms. Adding more number of H2 molecules would therefore decrease Ebinding still further. From Table 5 it can be seen that Ebinding of Li-decorated MOF-5 and Aldecorated MOF-5 lies between 18.00 kJ/mol and 11.70 kJ/mol and between 13.87 kJ/mol and 9.77 kJ/mol, respectively. These are much higher in magnitude than the usual pure MOF-H2 interaction strengths between 2.26e5.2 kJ/mol [6-8] mentioned earlier. Recalling the hydrogen binding energy range of ~24-34 kJ/mol estimated by Lochan and HeadGordan [9] at room temperature, we do not expect more than three H2 molecules per metal atom to remain adsorbed to Liand Al-decorated MOF-5s. Assuming that only sites near the metaldecorated BDCs are occupied and three H2 molecules per M atom implies 18 H2 molecules per formula unit of Zn4(BDC)3M6. This corresponds to a hydrogen uptake of 4.3 wt.% and 3.9 wt.% for Liand Al-decorated MOF-5, respectively. Thus, on basis of increased H2 binding energies and increased H2 storage capacities, we expect Liand Aldecorated MOF-5s to have potential for room temperature hydrogen storage applications.